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ABSTRACT 

The highly amphfied magnetic fields suggested by observations of some supernova remnant (SNR) 
shells are most likely an intrinsic part of efficient particle acceleration by shocks. This strong tur- 
bulence, which may result from cosmic ray driven instabilities, both resonant and non-resonant, in 
the shock precursor, is certain to play a critical role in self-consistent, nonlinear models of strong, 
cosmic ray modified shocks. Here we present a Monte Carlo model of nonlinear diffusive shock accel- 
eration (DSA) accounting for magnetic field amplification through resonant instabilities induced by 
accelerated particles, and including the effects of dissipation of turbulence upstream of a shock and 
the subsequent precursor plasma heating. Feedback effects between the plasma heating due to tur- 
bulence dissipation and particle injection are strong, adding to the nonlinear nature of efficient DSA. 
Describing the turbulence damping in a parameterized way, we reach two important results: first, for 
conditions typical of supernova remnant shocks, even a small amount of dissipated turbulence energy 
(~ 10%) is sufficient to significantly heat the precursor plasma, and second, the heating upstream 
of the shock leads to an increase in the injection of thermal particles at the subshock by a factor of 
several. In our results, the response of the non-linear shock structure to the boost in particle injection 
prevented the efficiency of particle acceleration and magnetic field amplification from increasing. We 
argue, however, that more advanced (possibly, non-resonant) models of turbulence generation and dis- 
sipation may lead to a scenario in which particle injection boost due to turbulence dissipation results 
in more efficient acceleration and even stronger amplified magnetic fields than without the dissipation. 
Subject headings: acceleration of particles - cosmic rays - supernova remnants - magnetic fields - 
turbulence - shock waves 



1. INTRODUCTION 

Observations of young supernova remnants (SNRs) 
suggest that strong collisionless shocks can simultane- 
ously place a large fraction of the s hock ram kinetic en- 
ergy i nto relativistic protons (e.g.. iBlandford fc Eichleij 
[19871 : iMalkov fc D^^ [20011 : IWarren et al.l 120051 ^ and 
amplify the ambient turbulen t magnetic field by large 
factors fe.g.JCowsik fc Sarkaj [T980: Reynolds & EUison, 
19921: iBamba et al.l 120031: Bcrczhko et al.l 120031: 



Vink k Laming! l2003t lUchivama fc AharonianI l2008f ). 
This coupling of diffusive shock acceleration (DSA) 
and magnetic field amplification (MFA) is critically 
important because the self-generated magnetic field 
largely determines the efficiency of DSA, the maximum 
particle energy a given shock can produce, and the 
synchrotron emission from radiating electrons. 

The generation and dissipation of strong MHD turbu- 
lence in collisionless, multi-fluid plasmas is a complex 
process. Different nonlinear approaches to the modeling 
of the large scale structure of a shock unde rgoing effi- 
cient c osmic ray acceleration (e .g., McKcnz ie fc Voelk 
19821: lAchterberg fc BlandfordI [198& iBell fc Lucek 
2OO1I : I Amato fc Blasil 120061 IVladimirov et alJ 120061 : 
I2OO6I : IZirakashvili et al.l l2008f l have 



predicted the presence of strong MHD turbulence in 
the shock precursor. However, an exact modeling of 
the shock structure in a turbulent medium, includ- 
ing nonthermal particle injection and acceleration, 
requires a nonperturbative, self-consistent description 
of a multi-component and multi-scale system includ- 
ing the strong MHD-turbulence dynamics. While a 
number of analytic models describing resonant and 
non-resonant amplification and damping of magnetic 
fluctuations have been proposed, these generally rely 
on the quasi-linear approximation that the fluctuations 
are small compared t o the background magnetic field, 
i.e., AB < Bq (e.g ., IWentzell [l97l iBvkov fc ToptvginI 
l2005t [KulsrudI [2005( 1. No consistent analytic description 
of magnetic turbulence generation with AB > Bq 
exists. For these reasons, numerical models with varying 
ranges of applicability have been proposed which offer 
a compro i nise between completeness and speed (e.g. . 
Belli [2OOI lAmato fc Blasl[200l [vladimirov et all [20061 : 



Zirakashvili et al.ll2008f) " 



Pelletier et al.l 



In principle, the problem can be solved completely with 
few assumptions and appr oximations with particle- in- 
cefl (PIC) simulatio ns fe.g.. [Bell[200l iSoitkovskvl [20081 : 
iNiemiec et al.l[2008[ ) or, in the assumption that electrons 
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are not dynamically important, by hyb rid models (e.g., 
IWinske fc Omidilll996HGiacalondl2004f ). However, mod- 
eling the nonlinear generation of relativistic particles and 
strong magnetic turbulence in collisionless shocks is com- 
putationally challenging and PIC simulations will not 
be able to fully address this problem in nonrelativistic 
shocks for some years to come even though they can pro- 
vide critical information on the plasma processes pro- 
ducing injection that can be obtained in no other way. 
In Appendix |^ we outline the requirements that a PIC 
simulation must fulfill in order to tackle the problem of 
efficient DSA with non-linear MFA in SNR shocks. 

In the Monte Carlo approximation we use here, the 
plasma interactions are parameterized allowing us to 
study coupled nonlinear effects between the extended 
shock precursor and the gas subshock. In particular, 
we investigate the nonlinear effects caused by upstream 
plasma heating due to magnetic field dissipation. 

The importance of the dissipation of turbulence in the 
shock precursor can be illustrated by the following es- 
timate. Suppose that in a shock wave of speed uq, the 
turbulence is generated by the resonant cosmic ray (CR) 
streaming instability, so the energy density of the tur- 
bulence, Uw, evolves approximately as u i^dUw{x)/dx — 
VAdPcr{x)/dx (e.g.. iBellfc Lucekl l200ll ). where Pcr(a;) 
is the CR pressure at position x and va is the Alfven 
speed. Ignoring all non-linear effects, the turbulence en- 
ergy density at the shock positioned at a; = is Uu]{0) = 
PoUqvao ■ Pci{0)/{poul). The ratio £„ = Pcr(0)/(/OoMo) 
characterizes the efhciency of acceleration and is typi- 
cally assumed to be on the order of ten percent or more. 
In the above, p is the fluid density and the subscript "0" 
always indicates far upstream values. Suppose a fraction, 
aH, of this energy goes into heating of the thermal gas 
in the shock precursor so the energy density of the ther- 
mal plasma increases by AUh{0) — anUwiO) at x = 0. 
Comparing AUh{0) with the internal energy density of 
the far upstream plasma, eo, we find 



eo Ma 



(1) 



where Ms is the sonic, and Ma the Alfven, Mach number. 
If rjH is large, the thermal gas in the shock precursor is 
strongly heated and this influences the subshock strength 
and the particle injection efficiency. In a non-linearly 
modified shock, a change in the injection efficiency causes 
the whole shock structure to change. For typical super- 
nova remnant (SNR) parameters (e.g., uq ^ 3000 km 
s-\ To ~ lO'^ K, no ~ 0.3 cm^^, and Bq - 3 ^G), the 
ratio M"^ /Ma ~ 250, and values of an as low as a few 
to ten percent may be important. 

Because existing analytical descriptions of MHD wave 
damping rely on the quasi-linear approximation AS ^ 
i?o, which is inapplicable for strong turbulence, and be- 
cause an exact description of this process in the frame- 
work of non-linear DSA is currently impossible (see Ap- 
pendix 1^ , we propose a parameterization of the turbu- 
lence damping rate. In doing this, we are pursuing two 
goals. First, we make some predictions connecting cos- 
mic ray spectra, turbulent magnetic fields and plasma 
temperatures, which, in principle, can be tested against 
high resolution X-ray observations in order to estimate 
the heating of the thermal gas by turbulence dissipation. 
And second, once heating is included in our simulation 



in a parameterized fashion, we will be ready to imple- 
ment more realistic models of turbulence generation and 
dissipation as they are developed. 

Our Monte Carlo simulation can be briefly summarized 
as follows (seelEl isMi et al.lll990aH Jones fc EllisCT3ll991l : 
IVladimirov et al1l2006L and references therein for more 
complete details). We describe particle transport in a 
plane shock by Bohm diffusion in a plasma flowing in 
the a:-direction with speed u{x). Particles move in small 
time steps as their local plasma frame momenta are 'scat- 
tered' at each step in a random walk process on a sphere 
in momentum space. Some shock heated thermal parti- 
cles are injected into the acceleration process when their 
history of random scatterings in the downstream region 
takes them back upstream. These particles gain energy 
and some continue to be accelerated in the first-order 
Fermi mechanism. This form of injection is generally 
called 'the rmal leakage' and w as first us ed in the cont ext 
of DSA in lEllisonet all (fT98l[ ) (see also lEllison|[l983 ). 

The magnetic field determining the random walk prop- 
erties through the diffusion coefhcient is the 'seed' (in- 
terstellar) magnetic field after it has been amplified by 
a large factor by the CR streaming instability and com- 
pressed and advected downstream with the flow. The 
streaming instability in the non-linear regime (AB 3> 
Bq) is described by the traditional quasi-linear equations, 
where the instability driving term is the CR pressure gra- 
dient. These quasi-linear equations are extrapolated into 
the non-linear regime in a parameterized fashion for lack 
of a more complete analytic description. The magnetic 
turbulence generated by the instability is assumed to dis- 
sipate at a rate proportional to the turbulence generation 
rate, and the dissipated energy is pumped directly into 
the thermal particle pool. An iterative scheme is em- 
ployed to ensure the conservation of mass, momentum, 
and energy fluxes, thus producing a self-consistent solu- 
tion of a steady-state, plane shock, with particle injec- 
tion and acceleration coupled to the bulk plasma flow 
modification and to the magnetic field amplification and 
damping. 

Our results show that even a small rate of turbulence 
dissipation can significantly increase the precursor tem- 
perature and that this, in turn, can increase the rate of 
injection of thermal particles. The nonlinear feedback of 
these changes on the shock structure, however, tend to 
cancel so that the spectrum of high energy particles is 
only modestly affected. 

2. MODEL 

The Monte Carlo simulation used here contains all o f 
the elements of the code used in IVladimirov et all (|2006D . 
and previously by Ellison and co-workers, i.e., it iter- 
atively determines a self-consistent, steady-state shock 
solution with resonant magnetic field amplification. The 
present code, however, has been completely re-written 
and optimized for the problem of DSA with MFA in non- 
relativistic, plane shocks. Besides including dissipation, 
the new code has been written for parallel processing and 
can model acceleration, in a reasonable time, to PeV en- 
ergies^. The particle transport is briefly described in 

^ The ability of the Monte Carlo code to accelerate pa rticles 
to PeV energies was shown in lEUison &: Vladimirovl 1)20081 ) . The 
results we are interested in here do not require such high energies 
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Section 12. li the amplification and dissipation of turbu- 
lence are described in Section 12.21 the effects of turbu- 
lence damping on the thermal plasma are described in 
Section [2. 3i and the iterative procedure used to reach a 
self-consistent solution is given in Section [2.41 

2.1. Particle Transport and Injection in the Monte 
Carlo Code 

Given a flow speed profile u{x) and a diffusion coef- 
ficient D{x,p) (which are obtained as explained below), 
the Monte Carlo code generates a thermal distribution 
of particles far upstream (or close to the shock, as de- 
scribed in Section 12.31 and Appendix [B]) , and propagates 
them by small time steps, St, scattering their momenta 
with t he 'pitch-angle scatter ing' scheme described in de- 
tail in lEUison et al.l (|1990af ). Since we assume strong 
turbulence {AB ^ Bq), the concept of 'pitch-angle' as 
an angle measured with respect to the average magnetic 
field direction loses it meaning. Instead, we define the 
'pitch-angle' as the angle that a particle's plasma frame 
momentum makes with the direction of the bulk flow. 
Each scattering is elastic and isotropic in the local plasma 
frame, which moves at speed u{x) with respect to the vis- 
cous subshock located at a; = 0, and the angle of scatter- 
ing is random, but its maximal value is determined by 6t 
and hy D{x,p) (see lEllison et al.lfl990al ). When particles 
cross the subshock, or move in any compressive how, the 
elastic scattering in the plasma frame makes them gain 
energy in the shock frame (see also Section fB.l.ip . 

Injection of particles into the acceleration process oc- 
curs in the Monte Carlo simulation when a formerly ther- 
mal particle first crosses the viscous subshock backwards, 
i.e., going against the flow. The number of particles that 
do this, and the energy they gain, are determined only 
by the random particle histories; no parameterization of 
the injection process is made other than the assumption 
of the diffusion coefficient value at various particle ener- 
gies.^ 

Since microscopic plasma modeling of particle injec- 
tion processes in non-relativistic subshocks with an ap- 
propriate 3-D PIC code is not yet available, all the mod- 
els dealing with particle acceleration by shocks must as- 
sume some prescription of the injection process. We fa- 
vor the present model for two reasons. First, our in- 
jection rate is set by the scattering prescription and 
doesn't require any additional assumptions. Second and 
even more important, our model was shown to agree 
well with spacecraft observations of the Earth's bow 
shock (Ellison & Moebius 1987; El lison etal.lll990bf ). in- 
terplanetary shocks (Baring et aLlll997), and with 1-D 



hybrid PIC simulations (Ellison et alj 



1991 . We are 



aware of altern ative models of pa rticle injection, such 
as that used bv iBlasi et al.l (|2005f l. in which the shock 
is assumed transparent only to particles exceeding the 
thermal gyroradius by a certai n factor. It i s pos sible to 
simulate the injection with the lBlasi et al.l (I2005D recipe 
in the Monte Carlo code, which was done bv lEllison et al.l 

and we use a smaller dyna mic range (see the note on dynamic range 
at the end of Section l3.2l l. 

^ As in previous implementations of our Monte Carlo model, 
the subshock is assumed to be transparent to all particles, includ- 
ing thermal ones, and possibly important plasma effects such as a 
cross-shock potential or large amplitude magnetic structures in the 
subshock layer are ignored. 



(|2005[ ). but it should be noted that an injection thresh- 
old may be i nconsistent w i th the b ow shock observations 
reported in lEllison et al.l ()1990bD . Our thermal injec- 
tion model is also simultaneously consistent with the bow 
shock helium and CNO observations with no additional 
parameters or assumptions. A comparative analysis of 
the different injection recipes is beyond the scope of our 
present study. 

A peculiarity of our approach is that in order to sepa- 
rate the CR particles from the thermal ones we use their 
history, and not their energy. By our definition, a ther- 
mal particle is one that we had introduced into the simu- 
lation upstream with a random thermal energy and that 
may have crossed the subshock going downstream, but 
has never crossed it back. Once a particle crosses the 
subshock (the coordinate a; = 0, to be more precise) in 
the upstream direction, it by our definition is injected 
and becomes a CR particle (see also Appendix [C|) . 

As particles are propagated and scattered, their con- 
tributions to fluxes of mass, momentum, and energy are 
calculated at select positions. We also calculate the pres- 
sure produced by the thermal particles Pthix) and the 
spectrum of pressure of the CR particles Pcy{x,p). (see 
details in Appendix [Bj) . The latter is related to the total 
CR pressure Pci{x) as 



(2) 



-Per (a;) = / Pciix,p)dp 



and Pci-{x,p) is then used to calculate the magnetic field 
amplification and dissipation, as described in Section 

We assume in this paper that the acceleration is size- 
limited, and model the finite size of the shock with a free 
escape boundary (FEB), located at position xfeb < 
far upstream of the shock. All CR particles crossing the 
boundary escape freely from the system. The maximal 
particle momentum, Pmaxj is thus determined when the 
upstream diffusion length becomes comparable to |a;FEB|- 
For a spherical SNR blast wave, |xfeb| is on the order 
of the radius of the remnant. 



2.2. Turbulence Amplification and Dissipation 

We assume that far upstream there exists a uniform 
magnetic field, Bq, parallel to the flow direction which 
is perturbed by transverse Alfvenic fluctuations with a 
power law energy spectrum. It is further assumed that 
these fluctuations produce Bohm diffusion for particles 
of all energies. Closer to the shock, where a popula- 
tion of accelerated particles drifting upstream is present, 
this seed turbulence is amplified by the CR streaming 
instability and additionally strengthened by the plasma 
compression. Once amplification begins, the spectrum 
of turbulence is no longer restricted to any particular 

^ We note that notatio n has changed slightly from that used 
in I Vladimirov et al.l II2006I ). The subscript 'tot' indicating integral 
quantities has been eliminated, and these quantities are now de- 
noted as Pct{x), L{x), etc. The quantities representing spectral 
densities of pressure, flux, etc., are denoted with the same letters, 
but a different set of arguments: Pct{x,p), L(x,k). To avoid ambi- 
guity, explicit definitions relating the spectral densities to the inte- 
gral densities are provided throughout the text. We also dropped 
the bar above the dissipation rate term L(x) for simplicity. 
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form^. Assuming that the turbulence is described by the 
quantities U-{x,k) and U+{x,k) {k is the wavenumber 
of turbulence harmonics, C/_ and [/+ are the spectra of 
energy density of structures propagating in the upstream 
and the downstream directions with respect to the ther- 
mal plasma, respectively), we model the evolution of the 
turbulence, as it is being advected with the plasma and 
amplified, with the following equations: 



E±[U] ^ {1 - aH)G±[U] + I±[U] 



(3) 



Here, for readability, we abbreviated as E the evolu- 
tion operator, as G the growth operator and as I the 
wave- wave interactions operator, acting on the spectrum 
of turbulence energy density U = {U-{x,k),U+{x,k)}. 
These quantities are defined as follows: 



E±[U]^{u±Vg)-^U± 



G±[U]-- 



I±[U]- 






dx\2 



-u±Vg 



Vgx 



dP„{x,p) 



dx 



dp 
dk 



'"go 



c/4 



(4) 

(5) 
(6) 



Here and throughout the paper, rgo = mpUoc/{eBo) - 
the gyroradius of a particle with a plasma frama speed 
equal to uq in the far upstream magnetic field Bq. The 
parameter an describes the turbulence dissipation rate, 
a nd for an = 0, equation ([3]) is exactly what was used 
in IVladimirov et all (|2006D , except there it was assumed 
that Vg <C u{x). Keeping Vq relative to u{x) results in 
a slightly greater amplification of magnetic field. In this 
system u = u{x) is the flow speed and Vg = Vg{x) is the 
parameter defining the turbulence growth rate and the 
wave speed^. 

For simplicity, we assume Vg{x) — Bq j ^J Att p{x) , 
in the present work, corresponding to /Aif = in 
IVladimirov et al.1 (|2006D . and emphasize that the use of 
equation ([3|) to describe the streaming instability when 
AS > Bq is only a parameterization. 

The growth operator G, which describes the turbu- 
lence amplification by the CR streaming instability, is 
proportional to the gradient of Pcr{x,p), the latter be- 
ing the spectrum of pressure of the CR particles driving 

* We assumed U± oc k^^ for the seed turbulence. The choice 
of the seed turbulence spectrum does not significantly affect our 
results for two reasons. First, the diffusion coefficient assumed here 
only depends on the total power in the turbulence and is insensitive 
to the shape of the spectrum, and second, the rapidly growing 
fluctuations due to the streaming instability quickly overpower the 
seed spectrum. Despite our using the Bohm diffusion coefficient, we 
still keep track of the turbulence spectrum in the simulation since 
this will be used in future work where the diffusion coefficient is 
determined self-consistently from the wave spectrum. 

^ As explained in IVladimirov et all I I2006I ). in the quasilin- 
ear case, AS <C Bq, the wave speed and the speed determin- 
ing turbulence growth rate are both equal to the Alfven speed, 
Vq(x) = Vj\ = Bq / yj A.-K p{x) . In the case of strong turbulence, 
AS > Bq, we hypothesise that the resonant streaming instabil- 
ity can still be described by equations JSjl with Vq being a free 
parameter ranging from Bq / yj4^Kp(x) to B^ff / ■^47rp(x). 



the instability. Pcr{x,p) is computed in the Monte Carlo 
simulation from the trajectories of particles, and the mo- 
mentum p at which it is taken in ([3]) is the momentum 
resonant with the turbulent structures with wavenumber 
k. The resonance condition assumed is 



eBo 



(7) 



The parameter an enters the equations of turbulence 
evolution ([3]) through the factor (1 — ctH), by which the 
term G describing the magnetic field amplification by 
the CR streaming is reduced (it is assumed that < 
an < 1). I3y writing this, we assumed simply that at all 
wavelengths only a fraction (1 — an) of the instability 
growth rate goes into the magnetic turbulence, and the 
remaining fraction an is lost in the dissipation process. 
The factor (1 — a n) in ^ can also be u nderstood in 
the following way: IVladimirov et all (|2006[ ) derived their 
equations (11) and (12) from equation (3) by assuming 
that the loss term L = {L here is L in their notation), 
but now we are assuming that 



L = -an ■ fa 



dx 



(8) 



For an = no dissipation occurs, and the CR stream- 
ing instability pumps the energy of the accelerated par- 
ticles into the magnetic turbulence, amplifying the ef- 
fective magnetic field most efficiently. For aH = 1 the 
additional turbulence energy produced by the instability 
is assumed to immediately dissipate, and the scattering 
in the shock is assumed to be provided only by the seed 
turbulence slightly increased by the plasma compression. 
The energy dissipation rate at all wavelengths is then 



00 
L{x) = / L{x, k) dk — an ■ Vg 



dPcrjx) 

dx 



(9) 



and we assume that all this energy goes directly into 
heating thermal particles. The modeling of the thermal 
plasma heating is covered in detail in Section 12.31 

When equation ^ is solved, the resulting U±{x, k) are 
used to calculate the amplified effective magnetic field 



Besix)= Att I U-{x,k)dk+ 



1/2 

47r / U+{x,k)dk 

the turbulence pressure 

BL(x) 



Pw{x) 

and the turbulence energy flux 



Stt 



(10) 



(11) 



Fw{x)=[-u^Vg] U^{x,k)dk + 



-u + Vg] / U+ix,k)dk 



(12) 
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see e quations (15), (18) and (19) in IVladimirov et all 
2006f )). which are then used in the derivation of a self- 
consistent solution, as discussed in Section \TM below. In 
the present paper we do not neglect Vg in the sum with 
u, which results in a slightly lower compression ratio rtot 
th an in the case 3/2n ± V g ~ 3/2u, that was adopted 
in IVladimirov et al.1 (|2006D . Note that the factor 3/2 
is valid for Alfven wave-like modes, which is implicitly 
assumed by our us ing the systeni of equ ations jS]). 

As mentioned by Caprioli ct al. (2008|) , effects from the 
transmission and reflection of Alfven waves at the sub- 
shock could be important, and accounting for these ef- 
fects may further lower the compression ratio rtot- We 
do not account for reflection and linear transformations 
of Alfven waves to other MHD modes (magnetosonic, en- 
tropy, etc.) at the subshock (see McKenzie and Westphal 
1970) in the present paper, because a correct description 
of these effects on the subshock in highly turbulent me- 
dia must contain simultaneously some other comparable 
effects. Indeed, MHD waves interacting with a shock pro- 
duce stochastic ripples in the shock surface and these rip- 
ples produce an effective broadening of the shock spatial 
structure de termined by the turbulence spectrum (see 
lBvkov|[l982f ). Moreover, as we argue in Appendix O 
accounting for suprathermal particles modifies the stan- 
dard Rankine-Hugoniot relations at the subshock and 
these effects could make it impossible to identify the sub- 
shock as a plane discontinuity, as assumed in all existing 
semi-analytic models. Clearly, these phenomena are im- 
portant and require further investigation but they are 
beyond the scope of our simplified description of MFA. 
Nevertheless, we believe our predictions regarding the 
turbulence dissipation in the shock precursor are quali- 
tatively correct. 

2.3. Heating of Thermal Plasma 
Repeating the deriv ation of equation (9) 



law suggest: 



iMcKenzie fc VoelM ()1982[ ). one obtains for a steady-state 
shock: 

"^ (Pti.p-'') = Hx) . (13) 



up' 



1- 



1 dx 



Here and elsewhere, p — p{x), u — u{x), and P^h = 
Pt\i{x). The quantity L{x) is the dissipation rate defined 
in ([H]), and the ratio of specific heats of an ideal nonrel- 
ativistic gas is 7 = 5/3. For L{x) = 0, equation (|13p 
reduces to the adiabatic heating law, Pth ^ p'^ and, for 
a non-zero L{x), it describes the heating of the thermal 
plasma in the shock precursor due to the dissipation of 
magnetic turbulence. The fiuid description of heating 
given by equation (fT^]) , while it doesn't include details of 
individual particle scattering, can be used in the Monte 
Carlo simulation to replace particle scattering and de- 
termine heating in the shock precursor. This merging of 
analytic and Monte Carlo techniques, or Analytic Pre- 
cursor Approximation (APA), is described in detail in 
Appendix [b1 

The APA only affects our treatment of thermal (i.e., 
not injected) particles in the precursor and involves two 
steps. The first is to introduce thermal particles into the 
simulation, not far upstream as we did before, but at 
some position xapa < close to the subshock and with 
a temperature equal to what Eq. (fT3l) and the ideal gas 



T{x) 



Pthix) 



ksnoiuo/uix))'' 



(14) 



where all quantities are taken at x — xapa and Pth is 
calculated from (|13[) (tiq and uq are the far upstream 
number density and shock speed, and fc^ is Boltzmann's 
constant). The second step is to calculate the thermal 
gas pressure throughout the precursor (at x < xapa) 
using (fTS]) instead of tracing particle motions. 

The main effects of turbulence dissipation in our 
model are: (i) a decrease in the value of the amplified 
field Bosix), which determines the diffusion coefficient, 
D(x,p); (ii) an increase in the temperature of particles 
just upstream of the subshock, which influences the in- 
jection of particles into the acceleration process, and (iii) 
an increase in the thermal particle pressure Pthix < 0), 
and a decrease in the turbulence pressure Pw{x), which 
enter the conservation equations described in Section ET^ 
Since all of these processes are coupled, a change in dis- 
sipation influences the overall structure of the shock. 

2.4. Closure of the Model 

Typically, we begin our simulation by propagating par- 
ticles, with a pre-set diffusion coefficient D(x,p), in an 
unmodified shock, where the fiow speed jumps discontin- 
uously from uq to U2, that is, u{x) — uq for a; < and 
u{x) — U2 for X > 0.^ This allows us to calculate the 
various ffuxes and other quantities, such as the CR pres- 
sure spectrum Pcr{x,p), at any position x. The latter is 
used to solve Eq. ([3]), which yields the turbulence spec- 
tra U±{x,k) and, subsequently, the amplified effective 
field Bae(x) and the pressure of the magnetic turbulence 
Pw(x). The spectrum Pcr(a;,p) also provides the turbu- 
lence dissipation rate ^ and the resulting pressure of 
the turbulence- heated gas (|13p . 

The equations for the conservation of mass and mo- 
mentum fiuxes are: 

p{x)u{x) = pqUq (15) 

(p and u are the mass density and the flow speed) and 

$p(x) = $po, (16) 

where $p(x) is the flux of the x-component of momen- 
tum in the x-direction including the contributions from 
particles and turbulence, and $po is the far upstream 
value of momentum flux, i.e., 

«'po - Poul + Ptho + Pwt) ■ (17) 

The quantity $p is defined as 

$p(x)= fp:,v,f{x,p)d^p + P^{x), (18) 



(here px and v^ are the a;-components of momentum and 
velocity of particles, and /(p) is their distribution func- 
tion, all measured in the shock frame), and in the sim- 
ulation it is calculated by summing the contributions of 

^ Everywhere in the text, unless otherwise noted, the subscript 
'0' indicates a far upstream value, '1' indicates a value just up- 
stream of the subshock, and '2' indicates a downstream value. For 
example, up = u (x = —00), mi = it(a:tr), and U2 = u{x > 0). See 
Appendix IB . 1 . ll for definition of xtr- 
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the particles crossing a grid location and adding the tur- 
bulence pressure P^ defined in (fTTj) . See details of this 
computation in Appendix [B] 

Initially, in our simulation, the shock doesn't have a 
self-consistent structure because we start with an un- 
modified shock and ^p{x) is overestimated at all loca- 
tions where accelerated particles are present. The next 
step is to use the calculated macroscopic quantities to 
find a new u{x) that reduces the mismatch between the 
local momentum flux and the far upstream value of it 
$P0 for a; < 0. We do this by calculating 



u'{x) = u{x) + s 



<^p{x) - $po 
PqUo 



(19) 



where u' {x) is the predicted flow speed for the next iter- 
ation, and s is a small positive number (typically around 
0.1), characterizing the pace of the iterative procedure. 
At this point we also refine our estimate for the particl e 
diffusion coefficient which, as in lVladimirov et al.l (|2006( ) . 
is assumed to be Bohm diffusion such that the particle 
mean free path is equal to its gyroradius in the effective, 
amplified field -Bcff(a;): 



D{x,p) = 



j\ 



vcp 



3eBes{x) 



(20) 



where Bce{x) is defined in pO)) . 

The predicted u{x) and D{x,p) are then used in a new 
iteration where particles are injected and propagated. 
The calculated CR pressure, momentum flux, etc. are 
then used to refine the guesses for u{x) and D{x,p) for 
the next iteration, and so on. This procedure is contin- 
ued until all quantities converge. 

In order to conserve momentum and energy, the com- 
pression ratio, rtot = uq/u2, must be determined self- 
consistently with the shock structure. To determine rtot 
we use the condition of the conservation of energy fiux 
given by: 



^e{x)+Qcsc{x) =$ 



EO, 



(21) 



where ^e{x) is the energy flux of particles and turbu- 
lence in the x-direction, Qesc is the energy flux of escap- 
ing particles at the FEB,'^ and the far upstream value of 
the energy flux is 



1 3 
*E0 = 2P0W0 



The quantity ^e{x) is defined as 



7 

H r-PthoWo + F^o- 

7- 1 



^e{x) 



Kv^f{x,p)(fp + F^{x), 



(22) 



(23) 



K being the kinetic energy of a particle with momentum 
p measured in the shock frame, and F^ is the energy flux 
of the turbulence defined in (fT2|) . The details of calcu- 
lating ^e{x) in the simulation are given in Appendix [Bl 
and the explanation of how ^^(x) is used in an itera- 
tive procedure converging to a consistent rtot is given in 
Appendix [Dj 

^ Particle escape at an upstream FEB also causes the mass and 
momentum flux es to change but these changes are negligible as 
long as uo < c l lEllisonlll985l ). 
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Fig. 1. — Dependence of magnetic field amplification and particle 
injection on the rate of magnetic turbulence dissipation in unmod- 
ified shocks. Results for shocks with compression ratios rtot=3.0, 
3.2, 3.4, 3.5 and 3.6 are represented by the thin solid, dotted, 
dashed, dash-dotted and thick solid lines, respectively. The x-axis 
variable is the turbulence dissipation rate an (constant throughout 
a shock), and the plotted quantities are the amplified downstream 
effective magnetic field -B(,ff2, the subshock sonic Mach number 
Msi, the fraction of simulation particles injected into the acceler- 
ation process /cr, and the ratio of amplified field to -Btrcnd down- 
stream from the shock. 



3. RESULTS 

3.1. Particle Injection in Unmodified Shocks (Subshock 

Modeling) 

In order to isolate the effects of plasma heating on 
particle injection, we first show results for unmodified 
shocks, i.e., u{x < 0) = uq and u{x > 0) = uo/^tot, with 
fixed rtot ■ In these models particle acceleration, magnetic 
field amplification and turbulence damping are included 
consistently with each other, but we do not obtain fully 
self-consistent solutions conserving momentum and en- 
ergy, since this requires the shock to be smoothed, while 
we intentionally fix u{x). 

In Fig. [1] we show results where the compression ratio 
is varied between rtot — 3 and 3.6 as indicated. In all 
models, uq = 3000 km s"\ Tq = 10^ K, no = 0.3 cm'^ 
and Bq = 3 /^G (the corresponding sonic and Alfvcn 
Mach numbers are M^o « Mao ~ 250). The FEB 
was set at 2:feb = —3 • 10^ rgo (our spatial scale unit 
rgo = TOpMoc/(ei?o)), and for each rtot we obtained re- 
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suits for different values of an between and 1. The 
values plotted in the top three panels of Fig. [T] are the 
amplified magnetic field downstream, BcS2, the Mach 
number right before the shock, Mgi (this is not equal 
to Mso because of the plasma heating due to turbulence 
dissipation) , and the fraction of thermal particles in the 
simulation that crossed the shock in the upstream direc- 
tion at least once (i.e., got injected), /„■ The bottom 
panel shows the ratio of the calculated downstream ef- 
fective magnetic field i?cff2 to trend values -Btrcnd(Q^i?); 
what is meant by "trend" is explained below. 

Looking at the curve for i3eff2 in the rtot = 3.0 model, 
one sees an easy to explain behavior: as the magnetic 
turbulence dissipation rate increases, the value of the 

amplified magnetic field decreases, going down to iJo'^tot 
(the upstream field compressed at the shock) for an = 1- 
Increasing an simply causes more energy to be removed 
from magnetic turbulence and put into thermal parti- 
cles, thus decreasing the value of B^m-^ In our model, 
the amount of dissipated turbulence energy scales lin- 
early with an- Therefore, the trend of i?cff2 changing 
with uh under the assumption that the total efficiency 
of the streaming instability is unchanged, but the en- 
ergy is channelled from magnetic turbulence to thermal 
particles, can be described as 

5trend(«i/) = (^0^?/') ^ + (1 " <^h) X 

~(Bo'i 

Off =0 V 



^cff2 



.3/4 V 
tot ) 



(24) 



where the first term on the right hand side is the Bq 
compressed at the shock (the compressi on factor r^/t i s 
explained by equation (13) in Vladimir ov et al.l I2006D . 
and the second term is proportional to the amount of 
the magnetic turbulence energy density generated by the 
instability for uh — 0, reduced by the factor (1 — uh)- 
Neglecting Bq, Eq. 



predicts Btrcnd oc \/l — ctfj- The 
comparison of B^fn with -Btrcnd(aff) from (j24p is shown 
in the bottom panel of Fig. [T] It is clear that the above 
trend matches the calculations very well for rtot =3.0. 
One also can see that the sonic Mach number in this 
simulation decreased from the upstream value of 250 to 
approximately 20, and that the fraction of injected par- 
ticles remained almost constant as an was varied. 

The curves for rtot = 3.2 demonstrate the same behav- 
ior, and the shape of the i?eff2 curve is similar to the one 
for rtot = 3.0, with the higher compression ratio produc- 
ing a higher value of the amplified magnetic field due to 
a greater number of particles getting injected. The cal- 
culated i3eff2 deviates from the trend i?trcnd(aff) more 
than in the rtot = 3.0 case, and this deviation marks the 
emergence of an effect that becomes more pronounced as 
rtot increases. 

The plots for rtot ^ 3.4 present a qualitatively dif- 
ferent behavior from those with rtot ^ 3.2. The down- 
stream magnetic field B^m does decrease with increasing 
an, but not as rapidly as in the previous two cases, and 
there is a switching point at an ~ 0.95 in the curves 
for Msi and /„■. The bottom panel of Fig. [T] shows a 
deviation of i?cff2 from the trend ((24|l by a large factor in 

* We use this result, with well understood behavior, to test the 
implementation of turbulence dissipation in our simulation. 



the rtot = 3.4 case. This effect becomes even more dra- 
matic for rtot = 3.5 and rtot = 3.6 where i3ef=f2j contrary 
to expectations, increases with an before an — > 1- The 
fact that the final energy in turbulence can increase as 
more energy is transferred from the turbulence to heat 
indicates the nonlinear behavior of the system and shows 
how sensitive the acceleration is to precursor heating. 

Indeed, if the upstream plasma is not heated suffi- 
ciently (as in the rtot < 3.4 cases), then the thermal 
particles first reaching the shock are cold (Ti ^^ Tq), and 
just upstream of the subshock the sonic Mach number 
Msi has a large value (much greater than 10 for any 
an)- As these thermal particles cross the shock, their 
momenta in the downstream plasma frame lie in a very 
narrow cone opening in the downstream direction, with 
the opening angle around 9 ~ AI~^ , making it equally 
difficult for most particles to turn around, cross the sub- 
shock backwards and get injected into the acceleration 
process. As long as 9 is small (or M^i is large enough), 
the number of injected particles is insensitive to the ex- 
act opening angle and fa stays relatively constant, as 
seen in the third panel of Fig. [T] for the rtot = 3.0 case. 

The injection rate increases with rtot and for rtot = 3.6, 
even with aH as low as 0.1, the CR-generated turbulence 
heats the thermal plasma through dissipation enough 
to lower Msi to '^ 10. Increasing an further lowers 
Msi even more, quickly increasing the probability that a 
downstream thermal particle will return upstream, thus 
boosting the injection rate (/„• rapidly goes up with an 
in the rtot > 3.4 cases). It turns out that this effect over- 
comes the reduction of i?cff2 due to damping, and BcS2 
starts increasing with an- 

For any rtot, at some high enough value of an (near 
0.90), the decrease in magnetic turbulence due to dissi- 
pation dominates the increase in injection and the mag- 
netic field drops. As an — > 1, and -Bcff2 becomes small 
enough, the efficiency of particle acceleration reduces suf- 
ficiently to cut down the total energy put into magnetic 
turbulence and the corresponding fraction of this energy 
dissipated into the thermal particles. This causes Mgi to 
turn up and /cr to turn down at about an = 0.95, as 
shown in the Fig. [T] 

It is worth mentioning that the observed increase of 
particle injection due to the precursor plasma heating 
is a consequence of the thermal leakage model of parti- 
cle injection adopted here. In this model, a downstream 
particle, thermal or otherwise, with plasma frame speed 
V > U2, has a pro bability to return upstream which in- 
creases with V (see lBelllll978l for a discussion of the prob- 
ability of returning particles). That is, we assume that 
the subshock is transparent to all particles with v > U2, 
but only some of these particles get injected, depend- 
ing on their random histories. Particles that don't get 
injected are converted far downstream out of the sys- 
tem. An al ternative model of injection (see, for example, 
iBlasi et aL 2005) is one where only particles with a gyro- 
radius greater than the shock thickness can get injected. 
The assumption with this threshold injection model is 
that the subshock thickness is comparable to the gyro- 
radius of a downstream thermal particle and only those 
particles with speeds v > ^z;th can be injected. Particles 
with V < S,vth are somehow blocked by the subshock. 
Here, fth is the downstream thermal speed and ^ is a 
free parameter, typically taken to be between 2 and 4. 
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Despite being similar, it can be expected that these 
injection models wi ll react differently to precursor heat- 
ing. Namely, in the lBlasi et al.1 (|2005f ) model the fraction 
of injected particles may be insensitive to the precursor 
heating if ^ is fixed^, because the same particle injection 
rate occurs regardless of the pre-subshock temperature 
Ti and downstream temperature T2. While both mod- 
els are highly simplified desc riptions of the com plex sub- 
shock (see, e.g. Malkov 1998i: iGiacalone fc Elli son 2000), 
they offer two scenarios for grasping a qualitatively cor- 
rect behavior of a shock where particle injection and ac- 
celeration are coupled to turbulence generation and flow 
modification. Hopefully, a clearer view of particle injec- 
tion by self-generated turbulence in a strongly magne- 
tized subshock will become available when relevant full 
particle PIC or hybrid simulations are performed. 

With the general trends observed here in mind, we now 
show how nonlinear effects modify the effect dissipation 
has on injection and MFA. 

3.2. Fully Nonlinear Model 

In this section we demonstrate the results of the fully 
nonlinear models, in which the flow structure, compres- 
sion ratio, magnetic turbulence, and particle distribution 
are all determined self-consistently, so that the fluxes of 
mass, momentum and energy are conserved across the 
shock. 

We use two sets of parameters, one with the far up- 
stream gas temperature To — W^ K and the far upstream 
particle density hq = 0.3 cm~^, typical of the cold inter- 
stellar medium (ISM), and one with Tq = 10^ K and 
no — 0.003 cm~^, typical of the hot ISM. In both cases 
we assumed the shock speed uq = 5000 km s~^, and the 
initial magnetic held Bq = 3 fiG (giving an equipartition 
of magnetic and thermal energy far upstream, nQksTf) « 
i?o/(87r)). The corresponding sonic and Alfven Mach 
numbers are Ms ~ Ma ~ 400 in both cases). The 
size of the shocks was limited by a FEB located at 
a^FEB = — lO^rgo ~ —3 • 10""* pc. For both cases, we ran 
seven simulations with different values of the dissipation 
rate an, namely an G {0; 0.1; 0.25; 0.5; 0.75; 0.9; 1.0}. 
Also, for the hot ISM {Tq = 10^ K) case we ran a sim- 
ulation neglecting the streaming instability effects, i.e., 
keeping the magnetic field constant throughout the shock 
and assuming that the precursor plasma is heated only 
by adiabatic compression (this model will be referred to 
as the 'no MFA case'). 



^ lAmato fc Blasil 120061) , who performed a brief analysis of the 
impact of the turbulence dissipation on the nonlinear shock struc- 
ture in a way similar to ours, did not report an influence of heating 
on the injection rate. 
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TABLE 1 

Summary of Non-linear Simulation in a Cold ISM 



an 


0.00 


0.10 


0.25 


0.50 


0.75 


0.95 


1.00 


''tot 


16.0 


16.2 


14.5 


14.6 


14.0 


13.2 


13.0 


'■sub 


2.95 


2.83 


2.75 


2.59 


2.50 


2.50 


2.51 


Beff2, MG 


345 


323 


284 


232 


158 


71 


21 


^trcndi MG 


345 


327 


299 


245 


174 


80 


21 


(r(x<o)>, 10-* K 


1.06 


4.3 


9.0 


17 


26 


37 


56 


Ti, 10* K 


3.3 


68 


160 


330 


490 


610 


650 


Ta, 10* K 


1400 


1500 


1600 


1600 


1800 


2000 


2200 


M,i 


44 


9.5 


6.3 


4.2 


3.5 


3.3 


3.2 


/cr, % 


1.0 


1.2 


1.6 


2.1 


2.6 


3.1 


3.2 


Pmax/n^pC 


500 


450 


400 


350 


250 


150 


80 


(7(a; < 0)) 


1.33 


1.33 


1.33 


1.33 


1.34 


1.34 


1.34 


72 


1.38 


1.38 


1.38 


1.39 


1.39 


1.40 


1.41 


3;tr/rg0 


-0.005 


-0.001 


-0.001 


-0.001 


-0.002 


-0.004 


-0.02 


I^APA/^'gO 


-0.04 


-0.06 


-0.06 


-0.09 


-0.23 


-0.57 


-2.1 



Note. — Here an, the fraction of dissipated energy, is the model 
input parameter, and the rest are results of the self-consistent simu- 
lation, as follows: rtot = uo/u2 is the total shock compression ratio, 
''sub = u\/u2 is the subshock compression, Sofia is the amplified effec- 
tive magnetic field downstream, -Btrcnd is the trend value calculated from 
Eq. (|24|l . {T{x < 0)) is the temperature of the precursor averaged over 
the volume from x = xfeb to a; = 0, Ti is the temperature at a:: = xtr, 
T2 is the volume-averaged temperature at a; > (all temperatures are 
calculated from the ideal gas law (|14|l ). M^i is the sonic Mach number 
at 2; = Xt-c, /cr is the fraction of injected thermal particles, pmax is the 
maximum particle momentum (i.e., the momentum at which /(p) starts 
falling off exponentially with p), (7(3; < 0)) is the value of the polytropic 
index of particle gas, calculated from particle pressure and internal en- 
ergy density as described in Appendix|Dl and averaged over volume from 
X = xfeb to a; = 0, 72 is the same quantity averaged over x > 0, Xtr 
is the point at which the subshock starts, as defined by Eq. (|B4p . and 
xapa is the point defined by Eq. HB11[I . at which thermal particles were 
introduced by the APA procedure as described in Appendix IB] 



TABLE 2 
Summary of Non-linear Simulation in a Hot ISM 



Oh 


0.00 


0.10 


0.25 


0.50 


0.75 


0.95 


1.00 


NoMFA 


rtot 


8.1 


8.2 


8.3 


8.0 


7.8 


7.4 


7.3 


13 


''sub 


2.92 


2.75 


2.55 


2.44 


2.22 


2.15 


2.12 


2.75 


Seff2, I^G 


62 


60 


55 


44 


32 


17 


14 


21 


■Btrond, ^J■G 


62 


59 


54 


45 


33 


19 


13 


- 


{T{x < 0)>, 10^ K 


1.04 


1.3 


1.7 


2.3 


3.1 


3.9 


4.2 


1.1 


Ti, 10^ K 


2.0 


6.0 


13 


23 


34 


42 


43 


2.7 


T2, 106 K 


53 


49 


47 


55 


62 


72 


75 


22 


M,i 


10.9 


5.8 


3.7 


2.6 


2.1 


1.9 


1.9 


4.7 


/cr, % 


1.2 


1.6 


2.5 


4.0 


6.4 


6.9 


6.4 


2.4 


Pmax /"IpC 


150 


120 


110 


100 


90 


70 


60 


80 


(7(a; < 0)) 


1.34 


1.34 


1.34 


1.34 


1.34 


1.34 


1.35 


1.34 


72 


1.43 


1.43 


1.43 


1.44 


1.45 


1.45 


1.45 


1.41 


aJtr/r-gO 


-0.04 


-0.02 


-0.02 


-0.02 


-0.03 


-0.07 


-0.05 


-0.02 


3;APA/''gO 


-0.1 


-0.1 


-0.2 


-0.2 


-0.4 


-0.9 


-1.4 


-0.1 



Note. — See the note for Table [T] for the explanation of listed quantities. 
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Referring to Tables [T] and [2l we summarize some of 
the results of these models. The effect of the turbulence 
dissipation into the thermal plasma is evident in the val- 
ues of the pre-subshock temperature Ti , the downstream 
temperature T2, and the volume-averaged precursor tem- 
perature {T{x < 0)) (the averaging takes place between 
X = xfeb and a; = 0). The temperatures were calcu- 
lated from the thermal particle pressure Pth{x) using the 
ideal gas law equation (fT4| . The value of Ti depends 
drastically on the level of the turbulence dissipation an, 
increasing from an = to an — 0.5 by a factor of 100 in 
the cold ISMcase, and by a factor of 11 in the hot ISM- 
case (less in the latter case, because the efficiency of the 
CR streaming instability in generating the magnetic tur- 
bulence is less for the smaller Mg and Ma for T = 10^ K). 
The values of the temperature as high as Ti are achieved 
upstream only near the subshock; the volume-averaged 
upstream temperature, {T{x < 0)), is significantly lower. 
The factors by which {T{x < 0)) increases in the above 
cases are 17 and 2.3, respectively. The value oi M^q/Mao 
is large in our models, and the estimate ([T]) predicts that 
even a small amount of dissipation is enough to raise the 
precursor temperature significantly. This is confirmed by 
our results: even an = 0.1 is enough to raise the pre- 
subshock temperature Ti approximately 20 times in the 
cold ISM (To = 10-* K) case. 

The downstream temperature, T2, varies less with 
changing an, because it is largely determined by the 
compression at the subshock, which is controlled by 
many factors as we discuss below. It is worth mention- 
ing the case without MFA reported in Table [51 Besides 
having a much larger compression factor than the shocks 
with MFA (rtot = 13 as opposed to rtot ^ 8), it has a 
much smaller downstream temperature (T2 = 2.2 • 10 K 
as opposed to T2 > 5.3-10^ K) These effects of dissipation 
on the precursor temperature may be observable. 

In Figure [2] we show results for /„, Mgi, i?cff2, and 
Strond which Can be compared to the results for unmod- 
ified shocks shown in Figure [T] For the modified shocks, 
the fraction of the thermal particles crossing the shock 
backwards for the first time, /„•, clearly increases by a 
large factor with an , which can be explained by the value 
of Msi dropping quickly below 10. One could expect that 
the amplified effective magnetic field -Bcff2 would behave 
similarly to the rtot = 3.5 case in Section l3.ll i.e. that 
BcS2 would not decrease or even would increase for larger 
an- Instead, i?cff2 behaves approximately according to 
the trend ([24]) , as the values of -Btrcnd from Tables [1] and 
[2| show and the bottom panel of Fig. [2| illustrates. The 
important point is that, even though precursor heating 
causes the injection efficiency to increase substantially, 
the efficiency of particle acceleration and magnetic tur- 
bulence generation is hardly changed. We base this as- 
sertion on the fact that i?cff2 remains close to -Btrcnd, 
which was derived under the assumption that changing 
an preserves the total energy generated by the instabil- 
ity, but re-distributes it between the turbulence and the 
thermal particles. The fact that the particle accelera- 
tion efiiciency is insensitive to an is directly seen in the 
results displayed in Figure [5| below. 

Considering how much the injection rate /d- increases 
with an, and how much the upstream temperature of the 
thermal plasma, Ti , is affected by the heating, it is some- 
what surprising that the trend of the amplified effective 
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Fig. 2. — Dependence of particle injection and magnetic field 
amplification on the rate of magnetic turbulence dissipation in the 
non-lincarly modified shock with u = 5000 km s~^ in the cold ISM 
and the hot ISM cases (the fully self-consistent models). The x-axis 
variable is the turbulence dissipation rate ajj (constant throughout 
a shock), and the plotted quantities are the amplified downstream 
effective magnetic field -B;,ff2, the subshock sonic Mach number 
Msi, the fraction of simulation particles injected into the acceler- 
ation process fcr , and the ratio of amplified field to -Btrend down- 
stream from the shock. 



field -Bcff2 is unaffected. The mechanism by which the 
shock adjusts to the changing heating and injection in 
order to preserve the MFA efhciency can be understood 
by looking at the trend of the total compression ratio 
rtot and the subshock compression ratio rsub in Tables [1] 
and [2 they both decrease significantly for higher an ■ 
The decrease in rsuh is easy to understand: with the tur- 
bulence dissipation operating in the precursor Mgi goes 
down, which lowers rgub- Additionally, decreasing P^i 
helps to reduce Tsub, and with a boost of the particle 
injection rate, the particles returning for the first time 
increase in number and build up some extra pressure 
just upstream of the shock, which causes the flow to slow 
down in that region, thus reducing the ratio rgub- The de- 
crease in rtot results from more complex processes. Here, 
the histories of 'adolescent' particles, i.e., those particles 
returning upstream for the first time or for the first few 
times, arc critical. Adolescent particles, while supcrthcr- 
mal, are still highly anisotropic in the shock frame and 
how they get accelerated in the smoothed precursor just 
upstream of the subshock determines the number and en- 
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Fig. 3. — Results of non-linear simulations in the cold ISM (left) 
and hot ISM (right) with different values of a^f . The solid, dashed 
and dotted lines correspond, respectively, to o^f = 0, 0.5 and 1.0. 
The plotted quantities are the bulk flow speed u{x), the effective 
amplified magnetic field B^f{{x) and the thermal gas temperature 
T{x). The shock is located at x = 0, and note the change from the 
logarithmic to the linear scale at x = —0.05 rgo- 



ergies of the 'mature' superthermal particles, i.e., those 
particles with enough energy to be nearly isotropic in 
the shock frame. The mature particles determine the 
CR pressure and precursor smoothing on larger scales. 

Further understanding of the shock adjustment to the 
changing dissipation can be gained by studying Figures 
[3] - [5l in which we plot the spatial structure and the 
momentum-dependent quantities of the shocks in the 
cold ISM and the hot ISM cases for an e {0; 0.5; 1}. 

Figure [3] shows an overlap in the curves for the flow 
speed u{x) in the an — and an — 0.5 models, and only 
close to the subshock u{x) falls off more rapidly towards 
the subshock in the oh — 0.5 case, resulting eventually 
in a lower rgub ■ This means that for the high energy par- 
ticles, which diffuse far upstream, the acceleration pro- 
cess will go on in about the same way with and without 
moderate turbulence dissipation (and the acceleration ef- 
ficiency will be preserved with changing uh)- For lower 
energy particles, however, there will be observable differ- 
ences in the energy spectrum. The an = 1-0 case has 
a significantly smoother precursor, which is not unusual, 
given the lower maximal energy of the accelerated par- 
ticles in this case (because of the magnetic field remain- 
ing low). The thermal gas temperatures T(x) plotted in 
the bottom panels of Figure [3] were calculated from the 
thermal pressure Pth{x) using (|14p and show that the 
temperature becomes high well in front of the subshock. 

In Figure 2] the subshock region for the hot ISM 
case is shown enlarged, and we can compare details of 
the models with (an = 0.5) and without dissipation 
(aH = 0.0). In the absence of turbulence dissipation, 
the thermal pressure Pth remains low upstream (see the 
middle panel), and the subshock transition is dominated 
by the magnetic pressure Pw For an = 0.5 (the bottom 
panel) the thermal pressure Pth just before the shock in- 
creases enough to become comparable with the magnetic 
pressure, but also the heating-boosted particle injection 
brings up the pressures of the 'adolescent' particles. As 
the plot shows, for an = 0.5 the pressures produced by 
the first and second time returning particles (Pi and P2) 
are not small compared to Pth and P^ just upstream 
of the shock, which contributes to the reduction of r^ub 
described above. However, the pressure of the 'mature' 
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Fig. 4. — Enlarged subshock region in the hot ISM case. The top 
panel shows the flow speed normalized to uq for the an = 0.0 and 
au = 0.5 models. The middle panel shows various constituents 
of pressure for the a^ = 0.0 run: the thermal pressure Ptjj, the 
pressures of 'adolescent' particles that crossed the shock 1, 2, 3 
and 4 times (Pi, P2, P3 and P4, respectively), the pressure of 
particles that have crossed 5 and more times P>5 and the magnetic 
turbulence pressure P^ , all of the above are normalized to the far 
upstream momentum flux <I>po. The bottom panel shows the same 
quantities for a[j = 0.5. 



particles, P>5, doesn't change much, which is a result 
of the non-linear response of the shock structure to the 
increased injection. 

The low energy parts of the particle distribution func- 
tions shown in Figure [5] are significantly different for 
models with and without dissipation in both the cold ISM 
and the hot ISM cases. The apparent widening of the 
thermal peak reflects the increase in the downstream gas 
temperature T2 . The differences extend from the thermal 
peak to mildly superthermal momenta 0.2 rripC, which 
shows an increased population of the 'adolescent' parti- 
cles with speeds up to u ~ 0.2c ~ 12mo when the turbu- 
lence dissipation operates. The high energy (p > 0.2TOpc) 
parts of the spectra for an = and an = 0.5 are similar 
(except for a lower Pmax due to a lower value of the am- 
plified field in the an = 0.5 case), confirming our asser- 
tion about the preservation of the particle acceleration 
efficiency. The increased population of the low-energy 
particles just above the thermal peak should influence 
the shock's X-ray emission. 

The characteristic concave curvature of the particle 
spectra above the thermal peak is clearly seen in the top 
panels of Figure [5l These shocks are strongly nonlinear 
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Fig. 5. — Results of non-linear simulations in the cold ISM(left) 
and hot ISM (right) with different values of aj{. Line styles as 
in Fig. |3] The plotted quantities are: the particle distribution 
function in the shock frame f{p) multiplied by p* (the normal- 
ization is such that J Airp'^ f(p')dp' = n, n being the number 
density in cm~^, and p' = p/(mpc)), the number of particles 
n(> p) with momentum greater than p (in cm"'') and the particle 
pressure (in dynes per cm^ ) per decade of normalized momentum 
dPp / dlogiQ p/ (mpc) . All quantities are calculated downstream at 
X = +6rgo ■ 



and, as the bottom panels in Figure [5] show, most of the 
pressure is in the highest energy particles. For these ex- 
amples, 60 to 80 percent of the downstream momentum 
flux is in CR particles. The number of particles produc- 
ing this pressure is small, however, and as the plots in 
the middle panels show, the fraction of particles above 
the thermal peak is on the order of 10^'^, and the fraction 
of particles above 1 GeV is around 10~^ in all cases. In 
addition to the pressure (and energy) in the distributions 
shown, a sizable fraction of shock ram kinetic energy flux 
escapes at the FEB. 

To summarize, for both the unmodified (Fig. [1]) and 
modified (Fig. ^ cases, Mgi drops and /„ grows as an 
increases. The surprising result is that i?off2 can increase 
in the unmodified shock as an goes up if rtot is large 
enough. This indicates that the boosted injection ef- 
ficiency (i.e., larger /„) outweighs the effects of field 
damping. This doesn't happen in the modified case (top 
panel of Fig. [2]) because of the nonlinear effects from the 
increased injection. From Fig. |3]we see that the boosted 
injection results in a smoother subshock and this makes 
it harder for low energy adolescent particles to gain en- 
ergy. Once particles reach a high enough momentum 
{p > 0.2mpc; see the top panel of Fig. [5]) they are able to 
diffuse far enough upstream where the boost in injection 
has a lesser effect. 

We must emphasise again that these results are very 
sensitive to the physics of particle injection at the sub- 
shock. It is difficult to predict how the nonlinear results 
would change if a different model of inje ction was used, 
but w e can refer the reader to the work of'Am ato fc Blasil 
(|2006n . who performed a similar calculation using the 
threshold injection model with a different diffusion coef- 
ficient. 

The free escape boundary in our simulation was rela- 
tively close to the shock (a^FEB = — lO'^r-go « — 3 x 10~^ 
pc), and the maximum accelerated particle energy was on 
the order of hundreds of GeV. These quantities, chosen 
to save computation time, are several orders of magni- 
tude short of typical SNR values. Nevertheless, we don't 



believe our results will be changed qualitatively if Pmax 
is increased. The reason is that even with our relatively 
low Pmax, the fraction of internal energy in relativistic 
particles is still large and the volume-averaged value of 
the polytropic index of the precursor plasma, (7(2; < 0)), 
as shown in Tables [T] and [21 is much closer to the 4/3 of a 
fully relativistic gas than to 5/3 for a nonrelativistic one. 
Increasing Pmax will not lower the polytropic index of the 
gas any further and, consequently, the plasma compres- 
sion and the subsequent acceleration efhciency will not 
change significantly (see Berezhko fc Ellisoniil99S ). 

4. CONCLUSIONS 

We have parameterized magnetic turbulence dissipa- 
tion as a fraction of turbulence energy generation and 
included this effect in our Monte Carlo model of strongly 
nonlinear shocks undergoing efficient DSA. The energy 
removed from the turbulence goes directly into the ther- 
mal particle population in the shock precursor. The 
Monte Carlo simulation self-consistently reacts to the 
changes in precursor heating and adjusts the injection 
of thermal particles into the DSA mechanism, as well as 
other nonlinear effects of DSA, accordingly. 

Our two most important results are, first, that even a 
small rate (~ 10%) of turbulence dissipation can dras- 
tically increase the precursor temperature, and second, 
that the precursor heating boosts particle injection into 
DSA by a large factor. The increase in particle injection 
modifies the low-energy part of the particle spectrum 
but, due to nonlinear feedback effects, does not signif- 
icantly change the overall efficiency or the high energy 
part of the spectrum. Both the precursor heating and 
modified spectral shape that occur with dissipation may 
have observable consequences. 

The fact that the shock back-reaction to the increased 
injection prevents the acceleration efficiency from chang- 
ing significantly is a clear consequence of the non-linear 
structure of the system. The boosted particle injection 
additionally smoothes the flow speed close to the sub- 
shock, which makes it harder for particles returning up- 
stream to gain energy. As a result, the population of the 
high energy particles is not much changed (except for the 
decrease in the maximum particle momentum due to the 
reduction in the effective amplified magnetic field from 
dissipation) and, because those particles carry the bulk of 
the CR pressure Per which drives the streaming instabil- 
ity, the amplification of the magnetic field is not strongly 
affected by the heating-boosted particle injection. 

The parameterization we use here is a simple one and 
a more advanced description of the turbulence damping 
may change our results. In our model the energy drained 
from the magnetic turbulence, at all wavelengths, is di- 
rectly 'pumped' into the thermal particles. Superthermal 
particles only gain extra energy due to heating because 
the thermal particles were more likely to return upstream 
and get accelerated. In a more advanced model of dis- 
sipation, where energy cascades from large-scale turbu- 
lence harmonics to the short-scale ones, the low energy 
CRs might gain energy directly from the dissipation. It 
is conceivable that cascading effects might increase the 
overall acceleration efficiency, the magnetic field ampli- 
fication, and increase the maximum particle energy a 
shock can produce. 

It is also possible that non-resonant turbulence insta- 
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bilities play an important role in magnetic field ampli- 
fication (e.g.. iPelletier et al.ll2006l ). This opens another 
possibility for the turbulence dissipation to produce an 
increase in the mae i ietic field amplification. For instance, 
iBykov fc Toptyginl (2005) proposed a mechanism for gen- 
erating long-wavelength perturbations of magnetic fields 
by low energy particles. If such a mechanism is respon- 
sible for generation of a significant fraction of the turbu- 
lence that confines the highest energy particles, then the 
increased particle injection due to the precursor heating 
may raise the maximum particle energy and, possibly, 
the value of the amplified magnetic field. 

While our model is for the most part phenomenologi- 
cal as far as particle transport, injection and acceleration. 



magnetic field generation, and dissipation are concerned, 
it allows us to investigate the coupled nonlinear effects 
in a shock undergoing efficient DSA and MFA. The more 
efficient DSA is, the more basic considerations of mo- 
mentum and energy conservation determine the shock 
structure and our model describes these effects fully. 
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APPENDIX 

A. COMMENTS ON PIC SIMULATIONS OF MFA 

There are two basic reasons why the problem of MFA in nonlinear diffusive shock acceleration (NL-DSA) is par- 
ticularly difficult for particle-in-cell (PIC) simulations. The first is that PIC simulations must be done fully in three 
dimensions to properly account for cross-field diffusion. As iJones et"al] (|1998D proved from first principles, PIC sim- 
ulations with one or more ignorable dimensions unphysically prevent particles from crossing magnetic field lines. In 
all but strictly parallel shock geometry,^" a condition which never occurs in strong turbulence, cross-field scattering is 
expected to contribute importantly to particle injection and must be fully accounted for if injection from the thermal 
background is to be modeled accurately. 

The second reason is that, in nonrelativistic shocks, NL-DSA spans large spatial, temporal, and momentum scales. 
The range of scales is more important than might be expected because DSA is intrinsically efficient and nonlinear 
effects tend to place a large fraction of the particle pressure in the highest energy particles (see Fig. [5]). The highest 
energy particles, with the largest diffusion lengths and longest acceleration times, feedback on the injection of the 
lowest energy particles with the shortest scales. The accelerated particles exchange their momentum and energy with 
the incoming thermal plasma through the magnetic fluctuations coupled to the flow. This results in the flow being 
decelerated and the plasma being heated. The structure of the shock, including the subshock where fresh particles are 
injected, depends critically on the highest energy particles in the system. 



A plasma simulation must resolve the electron skin depth, c/wpe, i.e., Lceii < c/tOpe-, where Upe 



[^-KTiee^ /rrieY'/'^ 



IS 



the electron plasma frequency and Lceii is the simulation cell size. Here, n,^ is the electron number density, rUe is the 
electron mass and c and e have their usual meanings. The simulation must also have a time step small compared to 
3, i.e., istop < ^pe ■ If "^6 wish to follow the acceleration of protons in DSA to the TeV energies present in SNRs 



pe ' 

we must have a simulation box that is as large as the upstream diffusion length of the highest energy protons, i.e., 
'«(Emax)/ito ^ ?'g(£'max)c/(3Mo), whcrc K is the diffusion coefficient, rg(£',„ax) is the gyroradius of a relativistic proton 
with the energy -Emax, uq is the shock speed, and we have assumed Bohm diffusion. The simulation must also be able 
to run for as long as the acceleration time of the highest energy protons, Tacc(^max) ^ E„^axc/{eBuQ). Here, B is some 
average magnetic field. The spatial condition gives 

At(-Emax)/^0 . 11 /i?max\ / Up VV ^ VV "e \ V 2/ / V^^ ( A\) 

(c/c,.) ^^^° l^^JUoOOkms-ij [-J^] iz—^) 1t^) > (Al) 



TcVy VlOOOkms^V VmG; Vcm-^y \1836J 
for the number of cells in one dimension. The factor / — mp/rrie is the proton to electron mass ratio. From the 
acceleration time condition, the required number of time steps is, 

6x1014 ( flE^ ) ( ^^L^V (£,]" f^^y/' (d^y^' . (A2) 
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LUpe \ TeVy VlOOOkms^V Kf^GJ Vcm-3/ V1836, 

Even with f ~ 1 these numbers are obviously far beyond any conceivable computing capabilities and they show that 
approximate methods are essential for studying NL-DSA. 

One approximation that is often used is a hybrid PIC simulation where the electrons are treated as a background 
fluid. To get the estimate of the requirements in this case we can take the minimum cell size as the thermal proton 
gyroradius, rgo = c^2mpEtii/{eB). Now, the number of cells, again in one dimension, is: 

7xl0M^^] f— ^^^^— t) 7^V'^' . (A3) 



K(£^max)/W0 



^kO 



IkeV 



The time step size must be istcp < ^cp^; where 
number of time steps to reach 1 TeV, 

'^acc v-^max 



TeVy V 1000 km s" 

p — eB jrnpC is the thermal proton gyrofrequency. This gives the 



1x10* 



Sn 



TeV 



Uq 



1000 km s" 



(A4) 



Parallel geometry is where the upstream magnetic field is parallel to the shock normal. 
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These combined spatial and temporal requirements, even for the most optimistic case of a hybrid simulation with an 
unrealistically large istop, are well beyond existing computing capabilities unless a maximum energy well below 1 TeV 
is used. 

Since the three-dimensional requirement is fundamental and relaxing it eliminates cross-field diffusion, restricting 
the energy range is the best way to make the problem assessable to hybrid PIC simulations. However, since producing 
relativistic particles from nonrelativistic ones is an essential part of the NL problem, the energy range must comfortably 
span nipC^ to be realistic. If E'max = lOGeV is used, with uq — 5000 km s^^, and Eth = lOMeV, equation (jA3p gives 
~ 1400 and equation (|A4p gives '^ 4x 10''. Now, the computation may be possible, even with the 3-D requirement, 
but the hybrid simulation can't fully investigate MFA since electron return currents are not modeled. The exact 
microscopic description of the system is not currently feasible. 

It's hard to make a comparison in run-time between PIC simulations and the Monte Carlo technique used here since 
we are not aware of any published results of 3-D PIC simulations of nonrelativistic shocks that follow particles from 
ful ly nonrelativi s tic to fully relativistic energies. A direct comparison of 1-D hybrid and Monte Carlo codes was given 
in lEUison et al.l (|1993f ) for energies consistent with the acceleration of diffuse ions at t he quasi-parallel Earth bow 
shock. Three-dimensional hybrid PIC results for nonrelativistic shocks were presented in iGiacalone &: EUisonl ([2000) 
and these were barely able to show injection and acceleration given the computational limits at that time. As for the 
Monte Carlo technique, all of the 16 nonlinear simulations presented in this paper were completed in approximately 4 
days on a parallel computing cluster, employing around 30 processors; enough statistical information was accumulated 
to restrict the uncertainty in the self-consistent value of rtot to about 5 percent. Increasing the dynamic range of 
the simulations to SNR-like energies (E'max ~ IPeV) would require 2-4 times more computation time. Thus, realistic 
Monte Carlo SNR models are possible with modest computing resouces. 

Despite these limitations, PIC simulations are the only way of self-consistently modeling the plasma physics of 
collisionless shocks. In particular, the injection of thermal particles in the large a mplitude waves and t i me varying 
struc ture of the subshock can only only be determined with PIC simulations (e.g., iNiemiec et al.l 120081 ISpitkovskvl 
|2008[) . Injection is one of the most important aspects of DSA and one where analytic and Monte Carlo techniques 
have large uncertainties. 

B. ANALYTIC PRECURSOR APPROXIMATION PROCEDURE 

Here we describe our Analytic Precursor Approximation (APA) where we introduce thermal particles into the MC 
simulation, not far upstream, but at some position xapa < close to the subshock. This procedure has two purposes. 
First, it saves computational time because we don't have to trace these particles along the extended shock precursor. 
Second, we use the APA to simulate the turbulence dissipation in the shock precursor: with this procedure we 
incorporate the analytic description of the heating of the thermal gas due to the turbulence dissipation into the Monte 
Carlo model of particle transport. 

In the next two subsections we describe how the momentum and energy fluxes in shock precursor are calculated in the 
absence of turbulence dissipation (part lB.l]) . and then explain how this scheme is changed when the Analyti c Pre cursor 
Approximation procedure is invoked to model the thermal plasma heating by turbulence dissipation (part IB.2[) . 



B.l. Precursor without Analytic Approximation 
B.1.1. Inherent Quasi- Adiabatic Heating and Subshock Definition in Monte Carlo Simulation 

It is worth pointing out that even if the turbulence dissipation is not included in our simulation, particles in the 
precursor will still be weakly heated. Just like an ideal collisionless gas put in a slowly shrinking volume will heat up 
adiabatically due to elastic collisions of particles with in-moving walls, the particles in our MC simulation traveling in 
a compressing shock precursor will heat up according to the adiabatic law Pth oc p''' ex u~'^ due to elastic scattering 
in the decelerating local plasma frame. This will be true as long as the particles have enough time to adjust to the 
changing flow speed. Quantitatively, the criterion for adiabatic heating is 

Tcoll < Tcomp, (Bl) 

where the collision time, TcoIi, is the ratio of the mean free path to the particle speed, (for nonrelativistic particles 
TcoU = fnpc/{eBcs)), and the compression time Tcomp is the temporal scale on which the speed of plasma that the 
particle is traveling in changes significantly, that is, by a value comparable to the particle speed v. Assuming u ^ w, 
Tcomp is approximately the ratio of the length on which the flow speed changes by v to the advection speed u: 
Tcomp = [v /\{du/ dx )\)/u . It also follows from the definition of Tcqh that the angular distribution of the particles 
remains isotropic if (jBip is true. The condition (jBip is therefore 



du 
dx 



« '-^^, (B2) 

urUpC 



restricting the flow speed gradient to small enough values, where adiabatic compression operates. 

If the magnitude of the flow speed gradient \du/dx\ is too large and (jBip doesn't hold, the increase of particle 
energies will be faster than adiabatic, and the angular distribution function will become non-isotropic, with more 
particles moving against the gradient than along it in the plasma frame. One obvious place where this occurs is the 
subshock. In our simulation, the subshock in the non-linear, self consistent solution gains a finite width, but the flow 
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speed jump at the subshock remains strong and rapid enough to efficiently accelerate particles. We define the point 
at which the transition from the adiabatic to the non-adiabatic regime occurs, xtr, by the following condition: 

Tcoll(a;tr) = ^Tcomp(2:tr), (B3) 



or, in terms of quantities available in the simulation, 

du 
dx 



1 veB^sjxtr) .-p,.^ 

6 u[Xtr)mpC 



where v = ■\/2fcBr(xtr)/"T-p, in which case |a;tr| is comparable to the local convective mean free path of a thermal 
particle. The factor 1/3 is chosen arbitrarily but our results are insensitive to it. Close to the shock, at xtr < x < 0, the 
flow speed drops so rapidly that the particles are heated in a non-adiabatic fashion. We can think of the non-adiabatic 
region as a subshock with a finite thickness \xf[\, and it is then reasonable to define the pre-subshock quantities 
(denoted with index '1') as the values at x = Xti-: wi = u{xtr), Tgub = u{xti-)/u{x > 0), etc. We note that Equation IB4I 
is only used to define the position of the subshock an is not used in any calculations. 

B.1.2. Direct Calculation of Momentum and Energy Fluxes 

The flux of momentum <^p[x) defined in (|18p is used in our simulation to calculate the smoothing of the precursor 
plasma fiow u{x), and the flux of energy $£;(a;) deflned in (|23p is used to calculate the compression ratio rtot consistent 
with the particle acceleration. The moments of the particle distribution function in psp and (|23D are the components 
of the stress tensor. If the plasma heating by turbulence dissipation is not modeled, and the Analytic Precursor 
Approximation procedure is not performed, then these moments are calculated in our simulation as described below. 
At select positions on the numerical grid spanning from far upstream to some position downstream from the subshock 
we sum the contributions of the particles that cross these positions to calculate the following: 

PxVxf{x,p)d^p = '^Pi^^Vi^^Wi, (B5) 

■' i 

j Kv.J{x,p)d''p = Y,K^Vi,^w^■ (B6) 

i 

Here the summation is taken over all particles crossing the position x at which the moments are calculated. The index 
i represents the considered particle, Px (v-x) is the a;-component of a particle's of momentum (velocity), K is the kinetic 
energy, all measured in the shock frame, and w is the weight of the particle defined as 



uo 



Vi,K 



t' <«^) 



In this definition the ratio |uo/''^i,x| is the weighting factor accounting for the fact that in our simulation particles 
crossing the position x at some angle to the flow do it less frequently than particles crossing parallel to the flow, no is 
the upstream number density of the plasma and Np is the number of simulation particles. 

If the particle distribution is isotropic in the local plasma frame moving at speed u{x) relative to the shock, then 
the quantities calculated in (jB5p and (jB6P can be expressed in the following way: 



PxVxf{x,p)d^p^p{x)u^{x)+Pp{x), (B8) 

Kvxf{x,p)d^p=-p{x)u^{x)+Wp{x)u{x) , (B9) 

where Pp{x) is the pressure and Wp{x) is the enthalpy of the particles. In the vicinity of the subshock, however, the 
isotropy assumption breaks down (see the discussion in Appendix IB.l.ll and the note on anisotropy in Appendix ICl). and 
the concept of isotropic pressure is not applicable. The fact that we directly calculate the moment of the distribution 
function in ([T8|) by evaluating the sum (jB5|l . instead of approximating the momentum flux with (jB8|l . ensures that 
we properly account for the effects of the anisotropy of particle distribution. This turns out to be important for 
self-consistently determining the flow speed u{x) near the subshock, which controls the subshock compression and the 
subsequent particle injection and acceleration efficiency. 

When plasma heating by turbulence dissipation is modeled, we replace the calculation (|18p with an analytic approx- 
imation assuming isotropic particle pressure, but we take special care to be certain that this approximation is only 
done far enough from the shock, where the isotropy approximation is applicable. This procedure is described below. 

B.2. Modeling Precursor Heating with the APA 

B.2.1. Particle Introduction Position 

The position at which the thermal particles are introduced must be close enough to the shock so that the analytic 
description of heating applies to most of the precursor extent, but far enough away from the non-adiabatic region 
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xtr < a; < 0, so that the analytic approximation remains valid where applied. In our simulation we chose the particle 
introduction position, xapa, so that the condition (jBip is only marginally valid at this position. Wc formalize it as 



du 
dx 



Tco\\{xapa) = MTcomp(a;APA), (BIO) 

which is equivalent to 

^ _^Z^eBeff(xAPA) ^ (Bll) 

u[XAPA)'mpC 



where we chose M — 0.1 and v = yj2k^jT^^/rrLp. At x < xapa we describe the thermal particle distribution function 
as a Maxwellian with the temperature defined by (fTH]) and (HU, and at a; > xapa we use the Monte Carlo simulation 
to describe the more complex particle dynamics. 

B.2.2. Momentum Space Distribution of Introduced Particles 

In order to include the effects of heating in the model, we must introduce thermal particles at xapa as if they were 
heated in the precursor, i.e., their temperature T(a::APA) must be determined by (fT3|) and (fT4|) . We therefore choose 
the magnitude of every particle's momentum p in the local plasma frame distributed according to Maxwell-Boltzmann 
distribution with probability density 

A / 1 \ 3/2 / 2 

4 / 1 \ 2 ..„ / P 



fip) = ^ U L. T^f V P^exp - " . (B12) 

0r \2TnkBT(xAPA) J \ 2mfcsT(xAPA) / 

The angular distribution of momenta of the introduced particles is a major issue of concern in doing simulation like 
ours because it determines the particle injection rate. We are replacing the dynamics of particles at x < a;APA with an 
analytical description, consequently we must distribute particles in p-space at a;APA the way they would be distributed 
having traveled from far upstream and reaching xapa for the first time. This is equivalent to calculating a p-space 
distribution of particles incident on a fully absorbing boundary at xapa after scattering in a non-uniform flow u{x). 
This is easy to do analytically if all particles have a plasma frame speed v less than the flow speed m(xapa) (because 
then all particles crossing position xapa do it for the first and the last time), and fairly complicated otherwise. We 
assume v < u{xapa) in further reasoning, which is justified by the fact that we find Ms{x apa) ^ 3 in most cases. 

As was stated earlier, we assume that the angular distribution of momenta of the introduced thermal particles is 
isotropic in the plasma frame. When these particles cross a position fixed in the shock frame, their fiux must be 
'weighted' to account for the fact that the number of particles arriving at xapa in a unit time is proportional to 
the cosine of the angle that their shock frame velocity Vgf makes with the x-axis. This can be done by assuming a 
probability density of /i = v^/v {v is the magnitude of the particle plasma frame velocity and v^ its x-component) as 

f{f^) = \{l + ^^l), (B13) 

where u — u(xapa)- It is normalized so that the probability Pr{jiQ < /i < /io + d/io) — fifJ'o)dfJ-o, and the functional 
form of (JBISP comes from the assumption that /(/i) oc fsf.x ~ u + /iu. 

B.2.3. Heating of the Upstream Plasma 

After the thermal particles are introduced at xapa and start to propagate in the shocked fiow, we have to calculate the 
momentum and energy fiuxes throughout the shock, for use in our iterative procedure. Because we didn't propagate the 
thermal particles at x < xapa, and because we must model the momentum flux redistribution between the turbulence 
and the thermal particles due to heating, we calculate the corresponding moments of particle distribution function the 
following way: 

I] Pi,xWi,xWi, if X > XAPA, 

p,vJ{x,V>)d'p^ { ;";)„2(^) ^ p^^(^) ^ J2 Pi,.V,,.W„ii X < XAPA, ^^^^^ 

ieCR 

J2 KiVi^^Wi, if X > XAPA, 

all i 

Kv,f{x,p)d^p= { Jp(x)u3(x) + ^_Pth(a:)u(x)+ (B15) 

+ J2 KiVi^^Wt,ii X < XAPA- 
leCR 

The summation for x > xapa is taken over all the particles crossing position x, and for x < xapa the summation index 
i e CR only includes the CR (i.e., injected) particles, while the contribution of the thermal particles is replaced by the 
analytic approximation of this contribution. The thermal pressure in this approximation is taken from the solution 
of p^ . For an = the equations (jB14p and (JBISP produce the same results (within intrinsic random deviations of 
the Monte Carlo code) as the calculation (|B5|1 and (jB6|l . 

In the region xapa < x < the heating due to L is ignored. This may, in principle, underestimate the heating of 
the upstream gas, but we find that this is a negligible effect. We prove it by running a simulation with another xapa, 
even farther away from the subshock, and making sure that the results are not affected significantly. 
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C. NOTE ON DIFFERENTIATION BETWEEN THERMAL AND CR PARTICLES 

As was mentioned earlier, we call a particle a thermal particle or a CR one based on its history: a CR particle 
is one that has gotten injected into the acceleration process by having crossed the shock from the downstream to 
the upstream region at least once^^. This criterion is used at two places in the model. First, it is used to calculate 
the spectrum of pressure driving the CR streaming instability Pcr{x,p): only the contribution of injected particles is 
included in Pcr{x,p). Second, when we calculate the thermal particle pressure Pthix) at a; < xapa using ^ and (|13p . 
and then introduce thermally distributed particles at a; = a;APA and continue the calculation of pressure for x > xapa 
from their trajectories as they elastically scatter in the flow, we implicitly assume that the dissipated energy of the 
turbulence goes directly into the thermal energy of the particles that have not been injected, i.e., thermal particles in 
our definition. 

If a description of particle-wave interactions in strong turbulence existed that explicitly described how particles of 
different energies participated differently in the instability generation and the turbulence dissipation, a criterion like 
ours would not be needed. However, such a description is not available and we believe that our way of separating thermal 
and superthermal (CR) particles for purposes of describing the instability growth and the turbulence dissipation grasps 
the essential non-linear effects in the shock structure. PIC simulations are, in principle, able to tackle this problem 
exactly but, as we mentioned earlier, they are extremely computationally expensive. 

We would like to point out an important consequence of our using the 'thermal leakage' model of particle injection 
into the acceleration process. Our simulation follows histories of charged particles from their 'childhood', when their 
speeds in the plasma frame are small compared to the bulk flow speed, to 'maturity', when they become relativistic. 
Unlike most semi-analytic descriptions of DSA, our model doesn't skip the 'adolescence' stage of particles, when after 
one or a few shock crossings the particles have speeds comparable to or slightly greater than the bulk flow speed. In the 
absence of these 'adolescent' particles, it is typically assumed that the jump in only the thermal particle pressure across 
the subshock determines the strength of the latter, and the superthermal part of the particle spectrum is continuous 
at the subshock and does not influence it. But the 'adolescent' particles that the Monte Carlo model does describe 
are not energetic enough to be insensitive to the subshock, and at the same time they have a strong anisotropy in the 
plasma frame, and therefore do not obey the Rankine-Hugoniot relations. This modifies the conservation laws across 
the subshock, because in the total kinetic pressure Pp{x) = Pth(a;) -I- -Per (2;) the term Pcr(a^) is not continuous at the 
subshock due to the contribution of the intermediate energy particles that it contains. 

D. COMPRESSION RATIO, TURBULENCE AND ESCAPING PARTICLES 

Equation (10) in Ellison e t al.l ()1990bl ) relates the fraction of energy flux carried away by escaping particles gesc to the 
total shock compression ratio, rtot- It assumes no magnetic field amplification and a polytropic index of downstream 
gas equal to 5/3 (in other words, neglects the effect of relativistic particles on the overall compressibility of the gas). In 
order to search for a rtot consistent with the shock structure and with particle escape, we derive a similar relationship 
that would account for the presence of magnetic turbulence and for the contribution of the relativistic particles. 
The problem is complicated by having particles of intermediate (mildly relativistic) energies and by the value of the 
magnetic field dependent on the particle acceleration. 

Writing equations p5|) . (fTH]) and ([2T|) for a point downstream of the shock, su fficie ntly far from it, so that the 
distribution of particle momenta is isotropic, and the approximations (jBSp and (|B9[) are valid, and denoting the 
corresponding quantities by index '2', we get: 

P2U2=PaUo, (Dl) 

P2ul + Pp2 + Pw2 = Poul + PpO + PwO = *P0, (D2) 

-P2M2 + Wp2U2 + -Fw2 + Qcsc = -^PoWq + WpQUQ + F^q = "^EO- (D3) 

The particle gas enthalpy Wp is Wp — ep + Pp, and the internal energy ep of gas is proportional to the pressure Pp. 
Introducing the quantity 7 so that Cp = Pp/{'y — 1), one can write 

w„u = _ ^ ^ P„u (D4) 

7 — 1 

The value of 7 is averaged over the whole particle spectrum, and it ranges between 5/3 for a nonrelativistic and 4/3 for 
an ultra-relativistic gas. The local value of 7 can be easily calculated in our code from the particle distribu_tion, along 
with Pp and Cp, as 7 = 1 -I- Pp/ep. Similarly, one can define S — F^/ {uP^) and calculate a local value of 5 anywhere 
in the code in order to express 

F^^6-P^u. (D5) 

For Vg <C u and Alfvenic turbulence, one expects ^ « 3, [see eq. ([T^ ]. 

Substituting (|D4p and (|D5P into the above equatio ns an d introducing rtot — uo/u2, we can eliminate p2 using (JD1[) 
and Pp2 using (|D2p . which allows us to express from (|D3P the quantity Qcsc = Qosc/^eo as 

.e.c = l + ^^^^;^^^, (D6) 

^^ Although the shock gains a finite width in the Monte Carlo simulation, we define the backward shock crossing as moving from x > 
to X < 0. 
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where 



A.^, (D7) 

72-1 

B^^^ (l + ^P" + ^"" ~ ^-A + g^, (D8) 



270 Ppo ^ 2^oPwO 
7o - 1 Po"o PoUq 



c^i + /^'\ ;p'; +i^i^. (D9) 



Note that Po^o/^pO = loM^, where 70 = 7 = 5/3 due to the absence of CRs far upstream, and, because we assume 

a seed turbulence far upstream, that provides a Bohm regime of scattering to all particles, one can write that far 

upstream AB w Bq, making PqUq/P^q = 2M^. 

The quantity Qcsc is readily available in the simulation after the end of any iteration. Comparing it to the value 

predicted by (jD6|l . we evaluate the self-consistency of the solution and make the correction to rtot, if necessary, for 

further iterations. For making these corrections it is helpful to use in the simulation the inverse of (jD6p . the physically 

relevant branch of which is 

2A 

rtot = , (DIO) 

B-^B^- 4AC{1 - fee) 



It is important to emphasise here that an iterative procedure similar to (|19p is still required to find the compression 
ratio rtot of a non-linearly modified shock, because quantities (/esc, Pw2 and 72 depend on rtot, so (jD10|) only provides 
a practical way to perform the iterations. 

E. NOTE ON SUBSHOCK PROPERTIES IN THE PRESENCE OF MFA 

iKang et al.l (|2002[ ) showed, using a thermal leakage injection model similar to what is used here, that in a non- 
linearly modified shock with Bohm diffusion and a sonic Mach number Mso, the subshock sonic Mach number scales as 
Msi ~ 2.9M^Q^^ with the corresponding subshock compression ratio Tgub = 4/(1 + 3/M^i). Numerically, for A/go = 400 
it gives Msi « 6 .3 and rgub ~ 3.7. 

The model of iKang et al.l (|2002l ) does not include magnetic field amplification and, as we have shown, in our case 
the values of the subshock sonic Mach number and compression ratio are quite different. For instance, for our model 
with MFA, but without magnetic turbulence damping {an = 0) we have found that Mso = 426 results in Mgi ~ 44 
and rgub ~ 3. The value of Mgi in the absence of precursor plasma heating is determined by the high CR pressure 
Per (a;), and the value of rgub in the high Mach number shock with MFA is largely controlled by the pressure of the 
amplified magnetic turbulence Pw{x) rather than by the thermal pressure Pth(a;)- The situation is changed when the 
turbulence dissipation operates, because it dampens Pw(x) and incre ases Pth(a^)i which reduces M^i- 

Our getting such a high value of A/gi is important because, just like lKang et al.l (|2002f l. we have found in Section [3Tl 
that Msi ~ 10 is a 'breaking point' in the thermal leakage model of particle injection: the injection rate depends 
weakly on Msi when Msi ^ 10, but increases rapidly with decreasing Afgi if M^i < 10. 

As an illustration for the above discussion, we show in Figure |6] the various constituents of the momentum fiux 
across the shock with an = in the cold ISM (Tq = 10^ K) case. It's clear that upstream of the shock the dominant 
contributor to the momentum flux is the CR pressure Per (a;) and, therefore, the precursor compression is mainly 
determined by Pcr(a;). For this particular set of parameters, Pcr(a;) results in a decrease in the flow speed u{x) by 
a factor of ~ 5.4 from its upstream value uq to the pre-subshock value ui. The temperature in the precursor only 
increases adiabatically for an = and this results in an increase over the far upstream temperature Ti/Tq « 3.3, thus 
reducing the local sonic Mach number to Msi ~ 44. 

The subshock compression ratio rgub is determined by the change in the different constituents of the momentum 
flux across the subshock. Pcr(a;), although large, changes little across the subshock^^, and what determines rg^b is the 
change in Pth(a;) and Pw{x). The latter, as the plots in Figure [HI show, contributes significantly to the momentum 
flux. This is an important point because the values of Poff2 and Pw{x) depend on rgub and rtot in a non-linear way 
making the traditional Rankine-Hugoniot relations inapplicable for determining rgub- Relation (JDIOP can be used to 
iteratively calculate the resulting compression ratio rtot, and it results in rtot — 16 and rgub = 2.95, as Table [T] shows. 
As we can see, the effect of magnetic turbulence pressure tends to decrease rgub compared to the case Pw{x) <C Pth(a;) 
(as in the latter case one would expect rgub ~ 4 for A/gi w 44). 



^^ Pct{x) is discontinuous at x = in our simulation due to the contribution to it of the particles that have crossed the shock only a 
few times. The jump in Pct{x) in this case is small compared to the jump in thermal and magnetic pressures across the subshock, but can 
be significant for ajj > 0, as seen in Figure |4] 
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X' ""go 

Fig. 6. — Illustration of rnorncntum flux balance in the non-linear simulation with a[f = in the cold ISM case. The thin dashed line is 
the thermal pressure P^i^{x), the thick dashed line - the CR pressure Per (a;), the dotted line is the magnetic turbulence pressure Pio{x), the 
dash-dotted line is the dynamic pressure p{x)u'^{x), and the thick solid line is the sum of the four, the total momentum flux. All quantities 
are normalized to the far upstream value of the total momentum flux. 
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